/*
 * point.h
 * This file is part of Jamming
 *
 * Copyright (C) 2012 - Yi
 *
 * Jamming is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * Jamming is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with Jamming. If not, see <http://www.gnu.org/licenses/>.
 */
#include<cstring>
#include <iostream>
#include <stdarg.h>
#include <exception>

namespace geo
{
/*class point*/
//exception
class NegativeOrZeroDimension : public std::exception {};

template< size_t nDim, class _coorType=int>
class point {
        _coorType coor[nDim];
    public:
        point(){
            if(nDim<=0) throw NegativeOrZeroDimension();
            memset(coor,_coorType(),sizeof(coor));
        };
        point(_coorType x,...) {
            va_list args;
            va_start(args,x);
            coor[0]=x;
            for(size_t i=1;i<nDim;i++) coor[i]=va_arg(args,_coorType);
            va_end( args );
        };
        point(_coorType* x) {
            memcpy(coor,x,nDim*sizeof(_coorType));
        }        
        //TODO implement range initializer 

           //operators
        point<nDim,_coorType>& operator +=(const point<nDim,_coorType>& p) {
            for(size_t i=0;i<nDim;i++) coor[i]+=p.coor[i];
            return *this;
        }
        point<nDim,_coorType>& operator -=(const point<nDim,_coorType>& p) {
            for(size_t i=0;i<nDim;i++) coor[i]-=p.coor[i];
            return *this;
        }
        point<nDim,_coorType>& operator *=(const _coorType& v) {
            for(size_t i=0;i<nDim;i++) coor[i]*=v;
            return *this;
        }
        point<nDim,_coorType>& operator /=(const _coorType& v) {
            for(size_t i=0;i<nDim;i++) coor[i]/=v;
            return *this;
        }
        const point<nDim,_coorType> operator +(point<nDim,_coorType>& p2) {
            return point<nDim,_coorType> (*this)+=p2;
         }
        const point<nDim,_coorType> operator -(point<nDim,_coorType>& p2) {
            return point<nDim,_coorType> (*this)-=p2;
         }
        //friend operators
        template< size_t nDim_, class _coorType_>
        friend bool operator!= (point<nDim_,_coorType_>& p1,point<nDim_,_coorType_>& p2);
        template< size_t nDim_, class _coorType_>
        friend bool operator== (point<nDim_,_coorType_>& p1,point<nDim_,_coorType_>& p2);
        template< size_t nDim_, class _coorType_>
        friend point<nDim_,_coorType_> operator* (const _coorType_& s,const point<nDim_,_coorType_>& v);
        template< size_t nDim_, class _coorType_>
        friend point<nDim_,_coorType_> operator* (const point<nDim_,_coorType_>& v,const _coorType_& s);
        template< size_t nDim_, class _coorType_>
        friend point<nDim_,_coorType_> operator/ (const point<nDim_,_coorType_>& v,const _coorType_& s);
        template< size_t nDim_, class _coorType_>
        friend std::ostream& operator <<(std::ostream& out,const point<nDim_,_coorType_>& p);
        template< size_t nDim_, class _coorType_>
        friend std::istream& operator >>(std::istream& in,const point<nDim_,_coorType_>& p);
        template< size_t nDim_, class _coorType_>
        friend _coorType_ dist(const point<nDim_,_coorType_>& p1,const point<nDim_,_coorType_>& p2);
        template< size_t nDim_, class _coorType_>
        friend _coorType_ operator*(const point<nDim_,_coorType_>& p1,const point<nDim_,_coorType_>& p2);
        //TODO normal to at least n-1 n-dimension vectors.
        template< size_t nDim_, class _coorType_>
        friend point<nDim_, _coorType_> normal(const point<nDim_,_coorType_>* pFirst,const point<nDim_,_coorType_>* pEnd);
    #define __point__type point<nDim_,_coorType_>
    //try use specification for nDim=2,3.
    //try use array 
    template< size_t nDim_, class _coorType_>
    friend point<nDim_,_coorType_> normal(const point<nDim_,_coorType_> p1,...);
    bool leftof(const point<nDim,_coorType> p1,...) {
        return ( normal(p1)*(*this)>0); 
    }
    bool on(const point<nDim,_coorType> p1,...) {
        return ( normal(p1)*(*this)==0); 
    }
}; // end of point
// Operation on point 
template< size_t nDim_, class _coorType_>
point<nDim_,_coorType_> normal(const point<nDim_,_coorType_> p1,...) {
            va_list args;
            va_start(args,p1);
            if(nDim_==2) return point<nDim_,_coorType_>(-p1.coor[1],p1.coor[0]);
            if(nDim_==3) {
                point<nDim_,_coorType_> p2=va_arg(args,__point__type);
                return point<nDim_,_coorType_>(p1.coor[1]*p2.coor[2]-p1.coor[2]*p2.coor[1],\
                                              p1.coor[2]*p2.coor[0]-p1.coor[0]*p2.coor[2],\
                                              p1.coor[0]*p2.coor[1]-p1.coor[1]*p2.coor[0]\
                                             );
            }
            //for(size_t i=1;i<nDim;i++) coor[i]=va_arg(args,_coorType);
            va_end( args );
            return p1;
};
template< size_t nDim_, class _coorType_>
bool operator== (point<nDim_,_coorType_>& p1,point<nDim_,_coorType_>& p2) {
        bool y=true;
        for(size_t i=0;i<nDim_;i++) y&=(p1.coor[i]==p2.coor[i]);
        return y;
};
template< size_t nDim_, class _coorType_>
bool operator!= (point<nDim_,_coorType_>& p1,point<nDim_,_coorType_>& p2) {
        bool y=false;
        for(size_t i=0;i<nDim_;i++) y|=(p1.coor[i]!=p2.coor[i]);
        return y;
};
template< size_t nDim_, class _coorType_>
point<nDim_,_coorType_> operator* (const _coorType_& s,const point<nDim_,_coorType_>& v) {
        _coorType_ newcoor[nDim_];
        for(size_t i=0;i<nDim_;i++) newcoor[i]=v.coor[i]*s;
        return point<nDim_,_coorType_>(newcoor);
}
template< size_t nDim_, class _coorType_>
point<nDim_,_coorType_> operator* (const point<nDim_,_coorType_>& v,const _coorType_& s) {
        _coorType_ newcoor[nDim_];
        for(size_t i=0;i<nDim_;i++) newcoor[i]=v.coor[i]*s;
        return point<nDim_,_coorType_>(newcoor);
}
template< size_t nDim_, class _coorType_>
_coorType_ operator*(const point<nDim_,_coorType_>& p1,const point<nDim_,_coorType_>& p2) {
        _coorType_ product=p1.coor[0]*p2.coor[0];
        for(size_t i=1;i<nDim_;i++) product*=p1.coor[i]*p2.coor[i];
        return product;
}
template< size_t nDim_, class _coorType_>
point<nDim_,_coorType_> operator/ (const point<nDim_,_coorType_>& v,const _coorType_& s) {
        _coorType_ newcoor[nDim_];
        for(size_t i=0;i<nDim_;i++) newcoor[i]=v.coor[i]/s;
        return point<nDim_,_coorType_>(newcoor);
}
template< size_t nDim_, class _coorType_>
std::ostream& operator <<(std::ostream& out,const point<nDim_,_coorType_>& p) {
         if(nDim_==1) out<<p.coor[nDim_];
         else {
             out<<"(";
            for(size_t i=0;i<nDim_-1;i++) out<<p.coor[i]<<",";
            if(nDim_>1) out<<p.coor[nDim_-1]<<")";
         }
         return  out;
};
template< size_t nDim_, class _coorType_>
std::istream& operator >>(std::istream& in,const point<nDim_,_coorType_>& p) {
         for(size_t i=0;i<nDim_;i++) in>>p.coor[i];
         return  in;
};
template< size_t nDim_, class _coorType_>
_coorType_ dist(const point<nDim_,_coorType_>& p1,const point<nDim_,_coorType_>& p2) {
        _coorType_ sum=_coorType_();
        for(size_t i=0;i<nDim_;i++) sum+=(p1.coor[i]-p2.coor[i])*(p1.coor[i]-p2.coor[i]);
        return sum;
    };
template< size_t nDim_, class _coorType_>
point<nDim_, _coorType_> normal(const point<nDim_,_coorType_>* pFirst,const point<nDim_,_coorType_>* pEnd){
};

/**
*  line 
*/
template< size_t nDim, class _coorType=int>
class line {
    point<nDim,_coorType> begin,end;
    public:
    line(point<nDim,_coorType> begin,point<nDim,_coorType> end) : begin(begin),end(end) {};
    //set up dx/dt,dy/dt,...
    void equation(_coorType par,...) {
    }
    point<nDim,_coorType>& origin() {
        return begin;
    };
    point<nDim,_coorType> direction() {
        return end-begin;
    }
    //operators
    
    //friend operator 
    friend std::ostream& operator <<(std::ostream& out,const line<nDim,_coorType>& l) {
        out<<l.begin<<"->"<<l.end;
        return out;
    }
    friend _coorType crossproduct(const line<nDim,_coorType>& l1,const line<nDim,_coorType>& l2) {
        return crossproduct(l1.end-l1.begin,l2.end-l2.begin);
    };
}; // end of line 

/*
* Point Set : implement with kd-tree
*/
class point_set {};

/*
* point space 
*/    
template< size_t nDim, class _coorType=int>
class space {
};
///* Class line */
//class line
//{
//    private:
//        /* data */
//    point<nDim,_coorType>& _beginTer,_endTer;
//    public:
//        line (point<nDim,_coorType>& p1,point<nDim,_coorType>& p2) {
//            _beginTer=p1;
//            _endTer=p2;
//        };
//        equation(_coorType par,...) {
//            //TODO set the equation of x=a_x*t+x_0, ...
//            // 1<-P_0, 2<- length of t, 3<-dX/dt...
//            for(_coorType* ptr=&x+1,size_t i=0;i<2*nDim;ptr-=2,i+=2) coor[i]=*ptr;
//        }
//        length();
//        
//        virtual ~line ();
//        
//}; //end of line
////operators
//point intersect(line l1,line l2);
//point parallel(line l1,line l2);


////Relationship between point & line
////TODO 
//// point is on line
//bool point::on(line l);
//// point is in the segment of line.
//point::in(line l); 
//point::above(line l);
//point::below(line l);
//// shortest distance to line;




///* Class polygon */
//class polygon
//{
//    private:
//        
//        
//    public:
//        polygon ();
//        virtual ~polygon ();
//}; //end of polygon
////operators of polygon


///* Class shape 
//*  a geometric shape defined by equation. 
//*/
//class shape
//{
//    private:
//        

//    public:
//        shape (arguments);
//        virtual ~shape ();
//}; //end of shape
////operators



}; /* Geo */



